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ABSTRACT 

Recent numerical and observational studies contain conflicting reports on the spectrum of magneto- 
hydrodynamic turbulence. In an attempt to clarify the issue we investigate anisotropic incompressible 
magnetohydrodynamic turbulence with a strong guide field Bq. We perform numerical simulations 
of the reduced MHD equations in a special setting that allows us to elucidate the transition between 
weak and strong turbulent regimes. Denote /cy, k± characteristic field-parallel and field-perpendicular 
wavenumbers of the fluctuations, and 6a the fluctuating field at the scale A ~ We find that 

when the critical balance condition, fc||i?o ~ k±bx, is satisfied, the turbulence is strong, and the en- 
ergy spectrum is E{k±) cx kj^^^ . As the fc|| width of the spectrum increases, the turbulence rapidly 

becomes weaker, and in the limit fc||i?o ^ k±b\, the spectrum approaches E{k±) cx kj^. The ob- 
served sensitivity of the spectrum to the balance of linear and nonlinear interactions may explain the 
confiicting numerical and observational findings where this balance condition is not well controlled. 
Subject headings: Magnetohydrodynamics, MHD, Reduced MHD, Turbulence 



1. INTRODUCTION 

In this Letter we investigate homogeneous and steadily 
driven incompressible magnetohydrodynamic (MHD) 
turbulence. In practical applications, such as fusion de- 
vices, solar wind, and interstellar medium, turbulence is 
driven by various large-scale instabilities, and turbulent 
energy is then spread over a broad range of spatial scales 
due to nonlinear interactions until small dissipative scales 
are reached where the energy is removed from the sys- 
tem. In the interval of scales between the injection and 
dissipation regions turbulence properties arc thought to 
be universal (e.g., Frisch 1995; Biskamp 2003). 

The MHD equations describing the evolution of mag- 
netic and velocity fluctuations b(x, t) and v(x, t) in the 
presence of a guide field Bq can be represented in the 
so-called Elsasser variables = v ± b: 

dtz^ T (Va ■ V)z± + (z^ • V)z± = -VP + f , (1) 

where Ya = Bo/\/47rp is the Alfven velocity, p is the 
fluid density, P is the pressure that is determined from 
the incompressibility condition, V • z='= = 0, f represents 
■ a large-scale forcing, and we omit the terms representing 
small viscosity and resistivity. The linear term on the 
left-hand side of equations H]), (V^ • V)z^, is responsi- 
ble for advection of the z"*" and z~ wave packets, with 
the Alfven velocity along the guide field. The nonlinear 
term, (z^ • V)z*, describes the interaction of turbulent 
fluctuations, and it is responsible for the energy trans- 
fer among different spatial scales. The nonlinear term is 
considered small if 

fc||Bo>fcj-&A, (2) 
where fc|| and k± are typical field-parallel and field- 
perpendicular wavenumbers of the fluctuations' spec- 
trum, and b\ (<C Bq) is the typical magnitude of fluc- 
tuations at the scale A ~ i/k±. This regime is referred 
to as "weak turbulence." 

The regime when the nonlinear term in not for- 
mally small will be called "strong turbulence." One 
can argue that in strong turbulence the following crit- 
ical balance condition should be maintained at all 



scales (|Goldreich fc SridhaiifTOQ^h : 

k\\Bo^kx^hx. (3) 

Indeed, during the characteristic time of nonlinear inter- 
action, ^ \/{k±b\), the fluctuations become corre- 
lated along the guide field up to a distance /|| ~ Vatn- 
This causality condition ensures the critical balance ([3]). 
Depending on the way turbulence is excited, it satisfies 
either condition ([2]) or ([3]) in a certain range of scales. 

Recent numerical simulations and analytic modeling 
suggest that in the case of strong turbulence ([3]), the 

field-perpen dicular en ergy s p ectrum is E(k\) o c fc , 
jM aron & Goldrcicfl [200H iMiiller fc Grappi^ [20051 : 
iBoldvrevI 12005. . 1200^ iMason et al. 1120061 l2007f) . How- 
ever, geophysical and astrophysical ob servations often 
exhibit somewhat stee per spectra (e.g., [Goldstein et al.l 
Il995t iBale. et all [20051 ). This raises the question of to 
what extent such systems can be described in the frame- 
work of MHD turbulence. 

In the present work we conduct direct numerical sim- 
ulations of reduced MHD equations, driven by a force 
with varying k\\ spectral width. This provides a unifying 
numerical setting allowing one to address the regimes 
of weak and strong turbulence in the same framework. 
We observe that when the critical balance ([3]) is satis- 
fied, the spectrum of strong MHD turbulence is close 
to —3/2. When the critical balance condition Q is 
even slightly broken, the spectrum steepens. As the 
weak turbulence condition ([2]) becomes better satisfied, 
the spectral exponent appro aches —2 in accord with the 
theory of weak turbulence l|Ng fc Bhattacharied [19961 : 
iGaltier et al,,2000 ) . The observed sensitivity of the spec- 
trum to the forcing details may explain conflicting re- 
sults of numerical and astrophysical observations, where 
the spectral properties of forcing are either not well con- 
trolled or not well known. 

Ac cording to the standard derivation (e.g., [Biskamnl 
[20031 ). the reduced MHD equations are valid in the re- 
gion k^_ ^ ; therefore, their applicability in the strong 
turbulence regime ([3|) is justified. Their applicability in 
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the weak turbulence regime however, requires an ex- 
planation, which we provide in the following sections. 

2. WEAK MHD TURBULENCE. 

When the condition ([2]) is satisfied, one may assume 
that turbulence consists of shear-Alfven and pseudo- 
Alfven waves, weakly interacting with each other. In 
the absence of nonlinear interaction, the waves would 
have random phases, and the Gaussian rule could be ap- 
plied to express their higher or der correlation funct ions 
through the second-order ones. 1 iGaltier et all (|2000l ) de- 
veloped a perturbative theory of weak MHD turbulence 
based on such random phase approximation. By expand- 
ing the MHD eq. ([T]) up to the second order in the non- 
linear interaction and using the Gaussian rule to split the 
fourth-order correlators, they derived a closed system of 
kinetic equations governing the wave energy spectra. 

These equations demonstrate that wave energy cas- 
cades in the Fourier space in the direction of large k±, 
and the universal spectrum of wave turbulence is estab- 
lished in the region k± 3> fc|| . In this limit the dynam- 
ics of the shear-Alfven waves decouple from the dynam- 
ics of the pseudo-Alfven waves, and the pseudo-Alfven 
waves are passively scattered by the shear-Alfven ones. 
The kinetic equation for the ener gy spectrum o f shear - 
Alfven waves, e(k,i), derived bv iGaltier et all (|2000( ). 
then reads: 



dte{k) = y"Mk,pqe(q)[e(p) - e(k)](5(g||)dk,pq 



(4) 



In this expression, the interaction kernel is Mk,pq = 
(7r/VA)(k_L X q_L)2(k_L • p±f /{k'lqlp'i), and we adopt 
the shorthand notation (ik,pq = ^(k — p — q)d^pd'^q. 
The phase- volume compensated energy spectrum is then 
calculated as E(k,t)dk^\ dk± — e(k,t)k±dk\\ dk±. The 
stationary solution of equation pi) was found ana lyti- 
cally and verified numerically in (|Galtier et al.|[2000l ). It 
has the general form £'(k) = where /(fc||) is 

an arbitrary function that is smooth at fc|| = 0. 

It should be noted, however, that the derivation of ([4]) 
based on the weak interaction approximation is not rig- 
orous. As follows from Eq. (jl]), only the q\\ — compo- 
nents of the energy spectrum e(q) are responsible for the 
energy transfer. However, if we apply Eq. ^ to these dy- 
namically important components themselves, that is, if 
we set fc|| = in 111 , we observe an inconsistency. Indeed, 
the perturbative approach implies that the linear fre- 
quencies of the waves are much larger than the frequency 
of their nonlinear interaction. The nonlinear interaction 
in Q remains nonzero as fc|| — > while the linear fre- 
quency of the corresponding Alfven modes, u>k = k\\VA -, 
vanishes. Therefore, as shown bv IGaltier et al.l (|2000[ ). 
the additional assumption of smoothness of the func- 
tion /(fc||) at /c|| = is crucial for deriving the spectrum 
E{k) oc kj^. 

A definitive numerical verification of such a spectrum 
seems therefore desirable. While numerical integration of 
a scattering model based on MHD equations expanded 
up to the second order in nonlinear interaction does 

^ Many papers contributed over the years to the develop- 
ment of fundamental ideas on MHD turbulence, see e.g., the re- 
views m idiskamp 2003; Ng ct al. 20Q3). The general methods 
of weak turbulenc e theory are reviewed in l lZakharov et al.1ll992l : 

EiiiireranilQoB). 



reproduc e the —2 exponent (jBhattacharjee &: Nell200"il : 
iNg et al.l [2003). this spectrum has not yet been con- 
firmed in direct numerical simulations of systems ([T]). 
The major problem faced by such simulations is to si- 
multaneously satisfy the two conditions, k± ^ fc|| and 
fc||-Bo 3> k±b\, which is hard to achieve with present-day 
computing power. In the next section we discuss a nu- 
merical setting in which the spectrum of weak turbulence 
can be verified. 

3. MODEL EQUATIONS. 

An important fact concernin g Eq. (U was emphasized 
bv IGaltier fc ChandranI (|2006f ). They noted that there 
exists a dynamical system that leads to exactly the same 
kinetic equation ([4]) in the weak turbulence regime ([2]), 
without any additional restrictions on fc_L and k\\ . To de- 
rive this system, we note that in the universal regime 
where Eq. ([4]) is applicable, the polarization vectors 
of the pseudo-Alfven modes are almost parallel to the 
guide field. One can therefore consider a system where 
such modes are eliminated for arbitrary k by restricting 
the initial MHD system to field-perpendicular fluctua- 
tions, z^: 

dti^ T (Va • V)z± + (zT . V)z± = -V^P + ^V2z± + f , 

Re 



= 0. 



(5) 



The fluctuating fields here have only two vector compo- 
nents, zi^ = {z^,5^,0}, but depend on all three spa- 
tial coordinates. Although system ([5]) is not presented in 
(jGaltier fc Chandranll2006D . their analysis of H]) is equiv- 
alent to solving such a system. Formally, syst em ^ is 
equivalent to the reduced M HD equations (e.g.. lBiskamd 
120031: IShebalin et all 119831 ). The principal difference is 
in the limits of validity: the reduced MHD model is ap- 
plicable only for fcj^ 3> fc|| , while we consider system ([5]) 
without any restrictions. 

Within the formalism of the weak turbulence theory 
both systems ([T|) and ^ lead to the same kinetic equa- 
tion (HI for the shear-Alfven turbulence. We thus con- 
clude that the wave energy spectrum obtained from the 
full MHD system ([T]) under the assumption k±^ ^ fcy 
should coincide with the spectrum obtained from the re- 
stricted system ([5]), where the condition k± ^ is not 
required. On the other hand, strong turbulence is ex- 
pected to develop when fc||-Bo ~ k±b\, which is precisely 
the domain in which reduced MHD provides a good ap- 
proximation of the full MHD model. This opens an effec- 
tive way for numerical investigation of both strong and 
weak MHD turbulence in the same framework, which is 
the goal of the present Letter. 

4. NUMERICAL METHOD AND RESULTS. 

We solve numerically the restricted MHD model (O 
using a fully dealiased pseudo-spectral technique in a pe- 
riodic box that is elongated along the guide field Bq with 
aspect ratio 27r : 27r : L^. The random force f has no 
component along z, it is solenoidal in the x — y plane 
and its Fourier coefficients outside the range 1 < fc^ < 2, 
[2t:/Lz) < < (2tt / Lz)nz are zero, where integer 
determines the width of the force spectrum in fcy . The 
Fourier coefficients inside that range are Gaussian ran- 
dom numbers with amplitude chosen so that the result- 
ing rms velocity fluctuations are of order unity. The in- 
dividual random values are refreshed independently at 
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Fig. 1. — The field-perpendicular energy spectrum E{k±) 
of MHD turbulence, calculated by direct numerical integration 
of l(5]l. The presented cases correspond, from bottom to top, to 
K = 1/6, 1/5, 2/5, 4/5, 8/5, 16/5. For clarity, the curves in the 
top panel are arbitrarily offset in the vertical direction. 

time intervals r = 0.1 L j_ / {2'KUrms) ■ The parameters 
and Lz control the degree to which condition ([2]) or ([3]) 
is satisfied at the forcing scale. Note that we do not 
drive the fc|| = modes but allow them to be gener- 
ated by nonlinear interactions. The Reynolds number 
is defined as Re = Urms{Li_/'2'K)/v, where L_l {— 27r) 
is the field-perpendicular box size, v is fluid viscosity, 
and Urms 1) is the rms value of velocity fluctua- 
tions. In our case magnetic resistivity and fluid vis- 
cosity are the same, v = rj. The system is evolved 
until a stationary state is reached, as determined by 
the time evolution of the total energy of the fluctu- 
ations. A typical run produces from 30 to 60 snap- 
shots. The field-perpendicular energy spectrum is ob- 
tained by averaging the angle-integrated Fourier spec- 
trum, Elk^) = 0.5(|v(kj_)j2)/c_L -t-G.5(|b(ki)|2)fci, over 
field-perpendicular planes in all snapshots. 

We performed a series of simulations for i?o = 5, = 
5L_L and Uz = 1, 2, 4, 8, 16. We used the resolution 
256^ mesh points in these simulations (and Re = 800), 
except for the case Uz = 16, where the resolution was 
512 X 512 X 256 (and Re±_ = 2000). We also performed 
a simulation with Bq = 5, Lz — 6L±, and riz = 1 
and the resolution 512 x 512 x 256 (and Re± = 2000). 
Fig. [T] shows the field-perpendicular energy spectra for 
each run. All the runs have different values of parameter 
K = {2Tr / Lz)nz that measures deviation from the critical 
balance ^ condition. 

We found that as the spectral width of the forcing 
along A;|| increases, higher and higher frequency modes 
of the velocity and magnetic fields are excited. For run 
K, = 1/6, all the forced modes have linear frequency 
oj = fc||_Bo ~ 1, which corresponds to a critically balanced 

forcing.^ In this case, the spectrum is E{k±) cx k^^^"^. 

^ It is important to keep in mind that the critical balance con- 
dition should be interpreted in a order-of-magnitude fashion, and 
that its validity is ultimately verified by the resulting spectrum of 
strong turbulence. 



Fig. 2. — Contour plots of the anisotropic energy spectra for 
the cases k = 1/6 (bottom), ft = 4/5 (center), and re = 16/5 (top). 
The colors represent energy in a log scale, normalized to the energy 
of the most dominant large-scale mode. 

This result is co nsistent with recent n umerical simulation 
of fuU MHD by iMason et al. I (|2007f) . since the reduced 
MHD system approximates the full MHD system when 
the critical balance condition is satisfied. As we increase 
the parameter k, we break the critical balance condition 
at the forcing scales. As a result, the spectrum monoton- 
ically steepens from —3/2 in the strong turbulence case 
to —2 in the weak turbulence shown in Fig. [TJ 

Fig. O shows isocontours of the full energy spectrum 
i?(fc||, fcj^) as a function of k± and for the three typi- 
cal cases K — 1/6, 4/5, 16/5. The bottom frame presents 
the energy distribution for the case k = 1/6, where the 
random force preserves the critical balance. As the cas- 
cade continues deeper into the inertial range, higher fre- 
quency modes co = k^^Bo are generated by virtue of non- 
linear interactions, just enough to maintain the critical 
balance condition at all scales, and establish a strong tur- 
bulence spectrum. As the frequency of the forced Alfven 
modes increases in the cases k — {2, . . . , 16)/5, the paral- 
lel cascade is slightly inhibited as the weaker interaction 
among the large scale Alfven modes dominates the en- 
ergy transfer to smaller scales, resulting in a steepening 
of the field-perpendicular energy spectrum. This can be 
seen in the middle and top frames in Fig. [21 where the 
distribution of energy becomes more and more elongated 
along k± rather than k^^ . 

5. DISCUSSION. 

Our numerical results demonstrate that if the energy 
spectrum has a limited extent, increasing the k\\ width 
of the forcing spectrum leads to the energy spectrum of 
weak turbulence c>c fc]]^. If however the width of the 
forcing is limited, but we can achieve arbitrarily high 
resolution in the k± direction, the interaction between 
Alfven modes will eventually become strong enough to 
satisfy critical balance and establish a strong turbulence 
spectrum. 

This is partly supported by the following derivation. It 
can be proved that turbulent fluctuations described by 
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system ([5]) satisfy the exact relation: 



(6) 



where Si^ = z='= (x + ) — (x) , and 5zf^ = • r /rj_ 
is the longitudinal component of . Averaging is taken 
over the statistical ensemble, or over time and spatial 
position X. In this formula, are the constant rates 
of z+ and z" energy dissipation. In the isotropic case, 
that is, without the g uide field, the relatio n analogous 
to (O was derived bv iPolitano fc Pouque"tl (fT998il . We 
now prove the following general inequality: 



any way; rather, the field-perpendicular energy spec- 
trum of strong MHD turbulence is flatter and closer 
to —3/2. This is cons is tent with recent re sults 
of iMiiller fc Grappinl (l200l : iMason et al.1 (|200l and 
also with hig h-resolution simulatio n s of isotropic MHD 
turbu lence by iHaugen et al.l (|2003D : iMininni fc Pouauel 
()2007( ). Astrophysical observations of the solar wind 
and of the interstellar medium reveal the presence of 
MHD turbulence, and find supp ort for both —5/3 and 
3/2 spectral exponents f e.g.. [Goldstein et al.l Il995t 



Goldstein fc R obert^ [l99l iBale. et all l2005t iBorovskvl 



2006; Podesta eTalJ (20061: ISmirnova et all I2006D . How 



{5z^{5i^ff<{{5i^f){{5i^T) 



(7) 



The first step is to use the Schwartz inequality, 
{5~zf{5i'^)Y < {iSz^)^){{Si^)^); the second step is to 
note that (Sz^)"^ < (dz^)'^, which completes the proof. 

Both sides in the expression ([7]) have finite limits as 
viscosity and resistivity go to zero. In the inertial in- 
terval, the correlation functions on the right-hand side 
of this expression have a power-law behavior, that is, 
{(Sz^)'^) oc and {{Sz^)'^) oc r'^ (we assume statistical 
symmetry between z"*" and z^). Since the left hand side 
of ^ scales as oc , and the inequality should hold for 
arbitrarily small r± , we obtain the exact result 

C2 + C4 < 2. (8) 

This inequality is useful for evaluation of these expo- 
nents from numerical simulations or experiments since 
the ratio of these exponents is wel l measured by th e 
method of extended self-similarity (jBenzi "eTaLlliggsl) . 
The inequality ^ then provides a boundary on the tur- 
bulence energy spectrum that is related to the second- 
order scaling exponent as E{k±) oc A;]|^~^^. In our case, 
the scaling exponent ^4 is usually close to 2^2 within 
small intermittenc y corrections , whic h can be checked 
numerically (e.g., IMiiller et al.l l2003f ). Inequahty ([8]) 
then implies that 3^2 < 2, and, therefore, the field- 
perpendicular energy spectrum cannot be essentially 

steeper than E(k±) oc kj_^^^ in the limit k± — > oo. 

Note, however, that in our numerical findings 
the spectral exponent —5/3 is not distinguished in 



ever, statistics of such data are often not good enough to 
distinguish between "—5/3" and "—3/2" with confidence. 
On the numerical side, simulations of MHD turbu- 
lence in the framewor k of reduced MHD \vere performed 
in m any works (e.g., iDmitruk et al.l [20031 : 1 Gomez et al.l 
I2005t iRapazzo et al.ll2Q07D : however, either the simula- 
tion domain was not anisotropic to ensure the critical 
balance condition ([3]), or the driving force was not spa- 
tially homogeneous, for example, applied at the bound- 
ary of the domain. 

Our results suggest that the interpretation of obser- 
vational and numerical results may be obscured if the 
fc|| and k± structure of the spectrum is either not well 
measured or not well controlled, in which case it is hard 
to deduce whether the field-parallel dynamics have been 
captured and whether the universal regime of MHD tur- 
bulence has been established. 
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